function pw = ODE(w, p ,par)
% Note: input p is a vector. w is a scalar.
investment = par.investment; investmentD = par.investmentD;  beta=par.beta;  alpha=par.alpha;  lambda=par.lambda;
AL=par.AL; AH=par.AH;  sigmaK=par.sigmaK; rho=par.rho; 
p_fitted = par.p_fitted;   p0=par.p0;  
print_level = par.print_level;

if(w>0 && w<1)
    % Solve jump size kappap
    psi = max( min(  (investment(p) + rho*p - AL)/(AH-AL),  1 ), w );   % xK >= 1 is imposed. A little bit regularization is very very important!
    xK = (psi/w);    yK=(1-psi)/(1-w);   
    
    p_post_jump = @(kappap)  ppval(p_fitted, w*(1-xK*kappap-alpha*max(beta*(xK-1), 0))./(1-kappap-w*alpha*max(beta*(xK-1), 0)) ) ;
    
    kappap = fsolve( @(kp) p - p_post_jump(kp) - kp , 0.01 ,  optimset('Display','off' ) );
%     for(  kappap_init = linspace( 0.001,  0.1, 4 ) )
%         kappap = fsolve( @(kp) p - p_post_jump(kp) - kp , kappap_init ,  optimset('Display','off' ) );
%         if( kappap>0  )
%             break;
%         end
%     end
    if(kappap<0 && print_level>=1)
        disp( ['kappap<0 when solving kappap at w=', num2str(w), ' and p=', num2str(p)] )
    end   
    
    sq = ( (AH-AL)./p -  lambda.*( kappap + alpha*beta ) ./ (  1-xK.*kappap  - alpha*beta.*(xK-1) ) + lambda.*kappap./(1-yK.*kappap)  ) ./ ( xK - yK )  ;
    
    if(sq<0 && print_level>=1)
         disp( ['sq<0 at w=', num2str(w), ' and p=', num2str(p)] )
    end 
    sigmap = sqrt( max(sq,0) ) - sigmaK;
    % Derivative pw
    pw =  max(0, sigmap./( w.*(1-w).*(xK-yK).*(sigmaK+sigmap) )) ;

elseif(w==0)
    disp( 'w=0 is evaluated in the ODE_noJump' );
    d0 = (AH-AL)/( investmentD(p0) + rho ) * (  (AH-AL)/( p0*sigmaK )^2 * p0 + 1  );
    pw = d0;     
else
    pw = 0;
end

